# Folders ####
directory<-Sys.getenv("US_Ineq_Repl")
CHwd<-file.path(directory, "Scripts/CH")
Datawd<-file.path(directory, "Processed")
Tabwd<-file.path(directory, "Results/Tables")

# Functions #### 
setwd(CHwd)
source("tablesCH.R")

# Code begins ####
setwd(Datawd)
load("CH0.RData")

mean0<-Mean
row.names(mean0)<-row.names(Mean)
Tab0<-table(Q1,Q2,Q3,Q4,Q5)

setwd(Datawd)
load("CH1.RData")
mean1<-Mean
row.names(mean1)<-row.names(Mean)
Tab1<-table(Q1,Q2,Q3,Q4,Q5)

setwd(Tabwd)
write.table(Tab0, 'CHPi0.csv', sep=',', row.names=T, col.names=T)
write.table(Tab1, 'CHPi1.csv', sep=',', row.names=T, col.names=T)

setwd(CHwd)
source("CHfun.R")

require(orthopolynom)
require(reldist)
require(matrixStats)
require(doParallel)

# Code begins ####
no_cores <- round(0.8*detectCores()) 
alpha<-0.05
y0<-1986
y1<-2015
D<-GetData(size=80000,boot=F,directory=Datawd,y0=y0,y1=y1)

W0<-D$Z0$lrwage3
W1<-D$Z1$lrwage3

Gini0<-Gini(W0)
Gini1<-Gini(W1)

Marg<-change(W0,W1,Boot=FALSE)

m<-1000
Rep<-1000
coe<-c('union','pubsect','manuf','nonwhite','female','partte','married','smsa',
       'exper','educ2','educ3','educ4','educ5')
start_time <- Sys.time()
marginals<-ParBootCI(R=Rep,nodes=no_cores,
                     BootW=BootW,Gini=Gini,Lorenz=Lorenz,BootZ=BootZ,
                     X_bar=X_bar,Gini_hat=Gini_hat,Ind_cov=Ind_cov,Ind_coef=Ind_coef,
                     mean0=mean0,mean1=mean1,
                     fun=change,W0=W0,W1=W1,m=m)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
(CImarg<-apply(marginals,1,quantile,probs=c(alpha/2,1-alpha/2)))
Aggr<-changeG(D$Z0,mean0,D$Z1,mean1,h=h,coe=coe,Boot=F)

start_time <- Sys.time()
aggregate<-ParBootCI(R=Rep,nodes=no_cores,
                     BootW=BootW,Gini=Gini,Lorenz=Lorenz,BootZ=BootZ,
                     X_bar=X_bar,Gini_hat=Gini_hat,Ind_cov=Ind_cov,Ind_coef=Ind_coef,
                     mean0=mean0,mean1=mean1,
                     fun=changeG,D$Z0,mean0,D$Z1,mean1,h=h,coe=coe,m=m)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')

eval(parse(text=paste0('Coe<-c("Cova","Coef",',paste(paste0('"',firstup(coe),'"'),collapse=','),',',
                       paste0(paste0('"',paste0(firstup(coe),'C'),'"'),collapse=','),')')))

CI<-NULL
for (cc in Coe) {
  eval(parse(text=paste0('CI',tolower(cc),'<-apply(aggregate[,colnames(aggregate)=="',cc,'"],1,quantile,probs=c(alpha/2,1-alpha/2))')))
  eval(parse(text=paste0('CI<-rbind(CI,CI',tolower(cc),')')))
}

T1<-100*cbind(Gini0,Gini1,Marg,Aggr)
colnames(T1)[1]<-"Gini0"
colnames(T1)[2]<-"Gini1"
colnames(T1)[3]<-"Change"
rownames(T1)<-c("25th","50th","75th","100th","Total")


T2<-100*t(rbind(CImarg,CI))
rownames(T2)<-c("25th","50th","75th","100th","Total")
T3<-apply(Aggr,2,'/',Marg)

Table<-matrix(rep("",3*nrow(T1)*ncol(T1)),nrow=3*nrow(T1),ncol=ncol(T1))
colnames(Table)<-colnames(T1)
rnames<-matrix(rep("",3*nrow(T1)),nrow=3*nrow(T1),ncol=1)
rnames[seq(1,3*nrow(T1),by=3)]<-rownames(T1)
rownames(Table)<-rnames

Table[seq(1,3*nrow(T1),by=3),]<-format(round(T1,digits=4),scientific=FALSE)
CIL<-matrix(rep("",nrow(T2)*ncol(T2)*0.5),nrow=nrow(T2),ncol=ncol(T2)*0.5)
colnames(CIL)<-colnames(Table)[3:length(colnames(Table))]
rownames(CIL)<-rownames(T2)

low<-T2[,colnames(T2)%in%"2.5%"]
ci_ni<-T2[,colnames(T2)%in%"2.5%"]<0
CIL[rownames(CIL)%in%rownames(ci_ni),][ci_ni]<-paste(format(round(T2[rownames(T2)%in%rownames(ci_ni),colnames(T2)%in%"2.5%"][ci_ni],digits=3),scientific=FALSE),";",sep="")
CIL[rownames(CIL)%in%rownames(ci_ni),][!ci_ni]<-paste(format(round(T2[rownames(T2)%in%rownames(ci_ni),colnames(T2)%in%"2.5%"][!ci_ni],digits=3),scientific=FALSE),";",sep="")

CIU<-matrix(rep("",nrow(T2)*ncol(T2)*0.5),nrow=nrow(T2),ncol=ncol(T2)*0.5)
colnames(CIU)<-colnames(Table)[3:length(colnames(Table))]
rownames(CIU)<-rownames(T2)

high<-T2[,colnames(T2)%in%"97.5%"]
ci_ni<-T2[,colnames(T2)%in%"97.5%"]<0
CIU[rownames(CIU)%in%rownames(ci_ni),][ci_ni]<-format(round(T2[rownames(T2)%in%rownames(ci_ni),colnames(T2)%in%"97.5%"][ci_ni],digits=3),scientific=FALSE)
CIU[rownames(CIU)%in%rownames(ci_ni),][!ci_ni]<-format(round(T2[rownames(T2)%in%rownames(ci_ni),colnames(T2)%in%"97.5%"][!ci_ni],digits=3),scientific=FALSE)

CI<-matrix(paste(CIL,CIU,sep=""),nrow=nrow(T2),ncol=ncol(T2)*0.5)
rownames(CI)<-rownames(CIL)
colnames(CI)<-colnames(CIL)

star<-matrix(rep("",nrow(low)*ncol(low)) ,nrow=nrow(low),ncol=ncol(low))
star[!(low<0 & high >0)]<-"*"
colnames(star)<-colnames(CI)
rownames(star)<-rownames(CI)

Table[seq(1,3*nrow(T1),by=3),colnames(Table)%in%colnames(star)]<-paste(Table[seq(1,3*nrow(T1),by=3),colnames(Table)%in%colnames(star)],star,sep="")
Table[seq(2,3*nrow(T1),by=3),colnames(Table)%in%colnames(CI)]<-CI
Table[seq(3,3*nrow(T1),by=3),4:length(colnames(Table))]<-format(round(T3,digits=3),scientific=FALSE)

Res<-matrix(rep("",3*nrow(T1)),nrow=3*nrow(T1),ncol=1)
Res[seq(1,3*nrow(T1),by=3),]<-format(round(T1[,"Change"]-T1[,"Cova"]-T1[,"Coef"],4),scientific=FALSE)
Res[seq(3,3*nrow(T1),by=3),]<-format(round((T1[,"Change"]-T1[,"Cova"]-T1[,"Coef"])/T1[,"Change"],3),scientific=FALSE)

Table<-cbind(Table[,1:5],Res,Table[,6:ncol(Table)])
colnames(Table)[colnames(Table)==""]<-"Res"

rownames(Table)[rownames(Table)==""]<-1:length(rownames(Table)[rownames(Table)==""])

setwd(Tabwd)
write.table(Table, 'CHtab.csv', sep=',', row.names=T, col.names=T)
setwd(Datawd)
rm(D)
save.image("CHtab.RData")
